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Abstract: The properties of the ground state of two-color QCD at non-zero baryon 
chemical potential /i present an interesting problem in strongly-interacting gauge theory; 
in particular the nature of the physically-relevant degrees of freedom in the superfluid 
phase in the post-onset regime /i > m, r /2 still needs clarification. In this study we present 
evidence for in-medium effects at high fi by studying the wavefunctions of mesonic 
and diquark states using orthodox lattice simulation techniques, made possible by the 
absence of a Sign Problem for the model with Nf — 2. Our results show that beyond 
onset the spatial extent of hadrons decreases as /j grows, and that the wavefunction 
profiles are consistent with the existence of a dynamically-gapped Fermi surface in this 
regime. 
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1 Introduction 


QC 2 D, the QCD-like gauge theory with gauge group SU(2), offers the simplest oppor¬ 
tunity for the controlled non-perturbative study of non-zero baryon charge density via 
Monte Carlo lattice simulation, unhindered by a Sign Problem. The reality of the quark 
determinant even with quark chemical potential /i^ 0 (and the consequent positivity of 
the path integral measure for number of flavors Nf even) arises from the pseudoreality 
of the fundamental representation of SU(2) [Ij. 

Several studies, both analytical and numerical, have helped elucidate the phase struc¬ 
ture and highlighted the differences with the physical case of SU(3). At low temperature 
T w 0, the ground state density n q of quarks rises from zero at a second-order onset 
transition at p — p 0 — permitting the existence of arbitrarily dilute baryonic mat¬ 
ter. For current quark masses sufficiently light to permit a separation of scales between 
m n and the next lightest hadron (either a 1 ~ p meson or 1 + diquark baryon), the regime 
p > p Q is described by chiral perturbation theory (yPT) [2j, and consists of a dilute 
gas of weakly-repelling Bose diquarks, which condense to form a superfluid ground state 
via the spontaneous breaking of the global U(l)# of baryon number. This regime has 
been identified in simulations with staggered lattice fermions mi- Higher densities 
have been probed in simulations with Wilson lattice fermions further from the chiral 
limit m- Over a wide range of p > p Q thermodynamic quantities such as n q , pressure 
p and energy density e |6] are all found to scale approximately with p as expected of 
a degenerate system of quarks, ie. one in which the available phase space is uniformly 
populated up to some characteristic momentum scale k F ~ /i, so that, eg. 
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in complete contrast to the behaviour predicted by yPT [2]. However, the system 
remains in the same confined state observed at p = 0, as indicated by the near-vanishing 
of the Polyakov line expectation ( L ). This state has some similarities with a confined 
chirally-symmetric ground state originally predicted in the context of large-W QCD [7]; 
for this reason we have called it a qua rkyonic regime. At the lowest temperatures 
studied, only for pa ~ 0(1) has color deconfinement signalled by (L) > 0 been observed; 
while the quark density n q is still much smaller than its saturation value 2NfN c /a 3 at 
this point, without simulations on a finer lattice we cannot yet be confident this high 
density deconfined phase persists in the continuum limit. Finally we note that the 
simulations suggest that the ground state is superfluid, as signalled by a non-vanishing 
diquark condensate, for all p > p Q ; however in the quarkyonic regime the scaling of the 
condensate ( qq(p )) is that expected from Cooper-pairing at a Fermi surface [5], rather 
than the characteristic BEC form found in y;PT ED- 

Beyond the fundamental requirements of determining the thermodynamic and sym¬ 
metry properties of the ground state, it is interesting to examine the nature of excitations. 
As well as offering continuity with the traditional concerns of lattice QCD at T — p — 0, 
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such questions bear on transport in the baryonic medium; answers to these questions 
in QCD would have the potential to inform, say, descriptions of neutron star spin down 
(via quantitative information on shear and bulk viscosities) and cooling (via a knowledge 
of which if any excitations remain gapless and hence capable of carrying energy away). 
There has been exploratory work in several directions. In [ 8 ] the hadron spectrum of 
QC 2 D was calculated as a function of p; beyond in the meson sector the usual ordering 
m-jr < m p is reversed, confirming earlier studies [9]. Above onset the lightest states are 
found in the 0 + and 1 + channels, with approximate degeneracy found between mesons 
and diquarks, as might be expected in a superfluid phase in which baryon number is 
no longer a good quantum number. The spectrum of heavy QQ quarkonium states also 
shows a non-trivial p-dependence | 10 ], possibly as a result of the formation of Qq states 
in the quarkyonic regime. In a recent study binding energies of multi-baryon “nuclei” 
formed from 0 + and 1 + bound states have been estimated [Ill- 

On a different tack, quark and gluon propagators have been calculated as func¬ 
tions of T and p in gauge-fixed configurations HUE]. The electric (longitudinal) gluon 
propagator in Landau gauge becomes strongly Debye-screened with increasing T and 
/ 1 , whereas the magnetic (transverse) gluon shows little sensitivity to T, and exhibits 
a mild enhancement in the quarkyonic regime before becoming suppressed at large p. 
Finally, the properties of topological excitations have been studied using a cooling pro¬ 
cedure to identify instantons ra- An enhancement of topological susceptibility \t is 
seen on entering the quarkyonic regime, which can be accommodated within the stan¬ 
dard perturbative description of Debye screening with the accompanying observation of 
a decrease in instanton scale size p(/i) oc /i~ 2 . Xt does fall very steeply, however, once 
( l ) > °. 

In this work we attempt to probe the interaction between quarks, and extract infor¬ 
mation on the spatial extent of hadrons, by calculating hadron correlation functions in 
which the qq or qq pair at the sink are spatially separated by a vector r m For a bound 
state H whose temporal decay in Euclidean space is governed by a simple exponential 
the spatial profile, determined numerically as a function of r. is proportional to 
the equal-time Bethe-Salpeter wavefunction 

T(r,r) = J d 3 T(O|'0(x, r)ip(x + r, t)\H). (2) 

The typical wavefunction profile for a bound state is gaussian, the width giving basic 
information about the size of the hadron. However, the correlators also yield interesting 
information even in the absence of a bound state, as explored in a study of the Z 2 
Gross-Neveu model with p 7 ^ 0 in 2+1 d |15]. Above onset, the wavefunction is no longer 
positive definite, but rather has an oscillatory structure with spatial frequency of order 
kp ~ p. These oscillations have a similar origin to the Friedel oscillations observed in 
the density-density correlations of electrons in metals (and thought to be responsible 
for the spin-glass behaviour of certain alloys), characteristic of a sharp, well-defined 
Fermi surface; the more primitive nature of the point-split hadron correlator makes 
it easier to measure in a numerical simulation, however. The observation of oscillatory 
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wavefunctions in [15] , with wavelength decreasing systematically with p, is one of several 
calculations leading to the identification of the Z 2 GN model as a Fermi liquid. 

A wavefunction study in QC 2 D has the potential to shed light on several outstand¬ 
ing issues in gauge theories at non-zero chemical potential, the most fundamental being 
whether it is indeed possible to identify a well-defined Fermi surface, since Fermi mo¬ 
mentum kp is not a gauge-invariant quantity. It may also help clarify the nature of 
the quarkyonic state, which roughly speaking may be thought of as a degenerate quark 
system in which only gauge-invariant excitations are permitted. Since two-quark inter¬ 
actions are the most relevant at a Fermi surface in the renormalisation group sense [16] , 
to what extent lessons learned with N c = 2 can be generalised to QCD remains to be 
seen. Nonetheless in principle the wavefunction should be a useful tool to chart the pas¬ 
sage from BEC to BCS realisations of superfluidity as p increases, which theoretically 
should take place for QC 2 D near enough the chiral limit. All these reasons motivate the 
current, exploratory study. 


2 Formulation 


In this section we explore the theoretical expectations for the wavefunction as a function 
of interquark separation r. We begin, following [IB], with the expression for the meson 
correlator C m (xo;r) with a local point source at the origin, and q and q separated by 
f at the sink. In anticipation of our later numerical results we choose the a priori 
arbitrary sign of p to yield the slowest decaying result in the positive Xq direction in 
diquark channels with non-zero baryon charge. Initially we assume free fields with 
quark mass m, and work at strictly zero temperature; the chemical potential p can then 
be understood as a Fermi energy for a system of degenerate quarks with Fermi energy 
Ep{ji) = p — \/kp + m 2 . The onset value at which the ground state contains a non-zero 
matter density is thus p Q = m: 

. r f l 4 r) r (j 4 a JP* -iqx -iq.r 

C m (z„,r) = X> j — j—T* ^ n + J ii _ n + m . (3) 


The Dirac matrix T = 1,75 for channels J p = 0 + , 0 . 

Performing an d the trace over Dirac indices, we obtain 


Cmi x o5 r) 


4 r r f 3 ^[±(p 0 + ip)(q 0 + ip) ± p 2 + m^e 1 ^ 0 ®Ao e ip-r 

(27 t ) 5 J q ° J P ° J P [(po + ip) 2 + P 2 + m 2 ][(q 0 + ip) 2 + p 2 + m 2 ] 

(4) 


where the ± sign denotes parity. First consider the case CF. The denominator of (JI]) 
has poles at p Q , q 0 = —ip ± iy/p 2 + rn 2 . Below onset, ie. for p < m, the integrals over 
Po and go ° n ly yield a non-vanishing result if the + root is chosen for po and the — root 
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for g 0 - We then use Cauchy’s theorem to obtain 


C m (x 0 -,r) = 


(2vr) 3 
2 

w? 


d 3 pe~ 2x ° -ip.f 


dpp 2 e - 2x °Vp 2 + m2 [ d6sm6e iprcos9 , 


(5) 


where in the last step r is taken to point towards the south pole. The angular integral 
then yields 

C m( x o',r) = ——twW dpph- 2x °VP 2+m2 Ji(pr), (6) 

(27r) 2 J o 

where Ji(z) = ( 2 / 7 tz)^ sin 2 is a Bessel function. Now expand (p 2 + m 2 ) in powers of 
p/m and use Laplace’s method of asymptotic expansion: 


C m (x 0 \r) » 


(2ir) 

jA 

(2ir) 


r -| e -2xom 


'0 


OO 2 

3 p xg 

dpp 2 e~ rn Jl(pr) 


m 

— I e 
x 0 


— 2mxo 


exp 


mr- 

4a; n 


(7) 


The corresponding expression for (7+ has an extra factor —p 2 /{jp 2 + m 2 ) under the final 
/7-integral, so that \C~\ > \C+\. The result (j7|) has engineering dimension 3 consistent 
with ([3]), is independent of p as expected below onset, and decays in euclidean time 
with mass 2m, also as expected. The wavefunction profile is a Gaussian whose width 
increases as y/x Q , as appropriate for two free particles gradually drifting apart. 

If instead we consider the corresponding cliquark correlator 


d 4 p f d A q 




(CF)t 


D ipx 


(27t) 4 + m 


(cry 


0 i qx p iq.r 


/i7o + m 


tr 


( 8 ) 


where the charge conjugation matrix satisfies C = — C 1 and C'yj/C 1 = — 7 ^, the result 
which emerges is 

Cb( x 0 , r ) = e 2fiX0 C*(x 0 , r). (9) 

Next consider the more interesting case above onset fi > m (the value of p 0 is 
governed by the mass per quark of the lightest baryon in the theory), where quark 
density n q = (g 7 og) > 0. For simplicity we set m = 0 so the only physical scale is // (for 
a sufficiently large Fermi surface this is always justified), and implying equality of + and 
— channels up to a sign. Now the po integral in Q can only receive a contribution from 
the pole at — i(p — \'p\) for \p\ > //, to yield 

2 POO 

Cm(x 0 -,r) = - —- 3 -r“3 / dpph~ 2pxo Ji(pr). (10) 

(27r)a 2 
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The only difference with (|6]) is the lower limit of the p-integral. On changing integration 
variable q = p — p and expressing the integrand as e~ 2x °^ q+ ^A(qg) where A can be 
expanded in powers of q using J' u (z) = J u -\(z) — (v/z)J u (z), we find 


C m (x 0 ;r) 


D -2x 0 p, 




(2t r 


dqe 


—2xo q 


'f 2 


i + 3 i 

2 fj, 


Ji (r/i) + qrJ[ (r/F) 


( 11 ) 


In fact, the leading behaviour as x 0 -+ oo is given by the term of O(q 0 ) 

1 


y P ~2px 0 


lim C m (x 0 : f) = 

Xq ->-00 


M 2e 


(2tt) 

2 


-Ji(r/i) (l + O(x 0 *)) 


(27r) 2 x 0 r 


Xo rs 

r: 2 '“" sin(r/i) (1 + t/ap)) . 


( 12 ) 


The expression (12) also has the correct engineering dimension. This time the state 
decays with energy 2/i, and arises from the promotion of a fermion from a negative 
energy state to the Fermi surface by injection of A E = 2p, A p = 0 (see Fig. [l]), 
leaving a negative energy hole which is re-interpreted as an anti-quark. The wavefunction 
envelope now decays algebraically (ie. oc r -1 ) rather than exponentially as in (Jr]) , and is 
modulated by oscillations of spatial frequency /i. Remarkably, the wavefunction can take 
negative values in this post-onset regime. The relation (J9]) continues to be respected; 
diquarks can be interpreted as gapless excitations (ie. A E = 0) formed from bound 
quarks with momenta ±p at the Fermi surface, as indicated in Fig. [l](c) . 

This is qualitatively similar to the phenomenon known in many-body theory as 
Friedel Oscillations , seen in the density-density correlation function in degenerate sys¬ 
tems. However, there are important differences; strictly Friedel oscillations have spatial 
frequency 2p, and decay with a different power [IS]. Both phenomena can be traced 
to the existence of a well-defined Fermi surface in the ground state above onset, lead¬ 
ing to a sharp lower limit on the integral in (10); the function C(xo,r) is much easier 
to measure using lattice techniques. The oscillations were observed in simulations of 
a 2+1 d IVy-flavor four-fermion model in [15]; here the applicability of free field theory 
could be justified by a calculation of the quantum corrections to the q - q interaction in 
the large-IV/ limit. In the static (ie. go = 0) limit this interaction is totally screened, ie. 
the effective Debye mass M(/i) —>■ oo. 

Since interactions between quarks clearly can’t be ignored in QC 2 D, a more am¬ 
bitious goal is to calculate C(x 0 -, r) in the presence of gauge interactions. As a first 
step, we consider one-gluon exchange between massless fermions using the Feynman- 
gauge propagator 5^F{q 2 ) (color factors are ignored). A simple possibility for F is 
the Debye-screened form (g 2 + M(/i) 2 ) -1 where in a gauge theory the expected relation 
is M/[x ~ 0(g) where g is the Yang-Mills coupling. The expression for C m is now a 
two-loop integral: 


C m = 9 2 ^2 tr 



F(q 2 


p Jp Jq 


.,ipx 


D —ip(x+r) 


A; 


lv 


i(i> + i ) - +7o ii> - +7o ip - F7o i(P + 


F7o 

(13) 
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Figure 1: Figure showing excitations of a system of massless degenerate fermions in the zero total 
momentum sector: (a) with meson quantum numbers requiring energy E = 2p at leading order (red); 
(b) with E <C /i following one-gluon exchange (green); and (c) with cliquark quantum numbers requiring 
zero energy (blue). Also shown in magenta is the modifying effect of a gap A > 0 opening at the Fermi 

surface. 


For simplicity’s sake, consider the static limit q 0 = 0 and T = 75 : 


C m (x 0 -,r) = J dp 0 J dp 0 j d 3 p J d^qe^-^e-^Fiq 2 ) x 

_ [(Po + ip){Po + ip) + p 2 } [(go + ip)(Po + ip) + (p + q) 2 ] _^ 

[(po + ip ) 2 + {p+ q) 2 } [(Po + ip ) 2 + P 2 ] [(po + ip ) 2 + P 2 ] [(Po + ip ) 2 + (p ■+ q ) 2 ]' 


The integrand of (14) has 4 poles po = —ip ± i\p + q\, po = —ip ± i\p\, and 4 
similar poles for p 0 , yielding in principle 16 cases to examine. However, the requirement 
for the Pq and po poles to be in opposite half-planes restricts us to only 6 non-vanishing 
possibilities. For either p 0 = —ip + i\p\ or p 0 = — ip + i\p + q] we can have po = —ip — i\p\ 
or po = —ip — i\p + q\; there will also be a contribution in a restricted corner of phase 
space from po = —ip+i\p+q\. It turns out that the contributions from the first two of the 
po poles listed decay in Euclidean time at least as fast as e~ 2fJ,x ° and e -A “° respectively, 
so we focus on the case p 0 = —ip + i\p\, po = —ip + i\p + q\- With this pole condition 
the integral becomes 


C m (x 0 -,r) 


4 g 2 

WY 



d 3 qe 


-(\P\~\p+q\)xo e ~ip-r 


F{q 2 


+ \p + q\) 


2 > 


(15) 


with the range of integration restricted to \p\ > / 1 , \p + q] < p , implying \p\~\pFq\ > 0 . 
There are regions of phase space, with |<f|/|pj m 0(1), where the integrand falls away 
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much more slowly than e~ px °. These regions are due to the excitation of a quark to a 
state lying above the Fermi surface leaving a hole in a positive energy state, with energy 
A E -C p in general, which is now permitted by the kinematics with A p = q ^ 0 (see 

I i- □ * 

Our tactic will be to perform f d 3 q first. The pole conditions dictate that the angle 
between q and — p is constrained to 6 q < O (p) = sin -1 and that the magnitude q is 
restricted to lie between pcos 6^(1 ± yA — cos 2 0/ cos 2 9 q ) (for 6 q = 0 this reduces to 
g 6 ip — p,p + p)). We arrive at 


x 


Ag 2 

W) 


p 2 dp 


dO p sin e p e iprcos6p 


' -©(p) 


dd q sin 6 q 


Qmax {jPfiq) 


dq- 


q2(p-yjp 2 +q 2 -2pq cos 6 q )x 0 


Iqm in (p,e q ) ( q 2 + M 2 )(2p 2 + q 2 - 2pq cos 6 q + 2p^/p 2 + q 2 - 2pqcos6 q ) 

V 


(2n 


_ 1 
'f 2 


dpp 2 Ji (pr)/(p, xq), 


(16) 


where the Debye-screened form of F(q 2 ) has been inserted. 

To examine things more closely, consider the limit 6 q = 0 so that p + q is either 
parallel or anti-parallel to p. The integral over q in (16) then becomes 


dq- 


q2 e ~qxo 


CP+P 


+ 


dq 


?-qx o 


q 2 + M 2 ' 


Jp _p {q 2 + M 2 )(2p — q)‘ 

After some rearrangement the integrals may be performed to yield 

e ~(p+p)xo e ~(v-p)xo(p _ nY 


(17) 


1 

x 0 


(p + p) 2 + M 2 


+ e 


-px 0 


x 0 + 


{p + p) 2 [(p - p) 2 + M 2 ] 


(1 + 0 (xq 1 ). (18) 


Reassuringly the factor multiplying e~ px ° cancels at this order in Xq this dependence 
could only have arisen via incomplete cancellation between the inner limits on the in¬ 
tegrals in dl7| ), but the behaviour of the integrand at the origin of p-space should be 
smooth. Leaving aside for now the slightly ill-defined integral over 6 q , we arrive at an 
approximate expression of the form 


()Tq , r) 


Ag 2 


R)(2tt) 5 
Ag 2 


p 2 dp 


(P P i ) x 0 (^ ^ 

dd sin 9 e ipr cos0p _ — _—_ 

P V ( p + M )2[(p_ M )2 + M 2] 


(2tt 
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i 2 £or 2 


dpp 2 Ji (pr 


ip — p) 2 e ^ ^ x ° 
ip + p) 2 ip - p) 2 + M 2 ’ 


(19) 


where terms decaying faster than e 2px ° and/or x 0 2 are neglected, and in the second 
line the integral over 9 p has been completed. 

















Finally, as in the one-loop case we substitute q — p — p and Taylor-expand the 


integrand in powers of q. The identity / 0 °° q n e qx °dq = n\/x q +1 is exploited to yield 


Cm(x Q ,f) 


2 g 2 J' 2 (H 1 

(27t) 8 (/ir)2 XqM 2 


(1 + O(x 0 )) 


Xn 


-2 


Ji (rp) 


(Xofi) 2 (r/j) 


( 20 ) 


The engineering dimension is wrong due to the approximation that has been carried out 
preceding eqn. There is no exponential decay in xq, resulting from the presence in 

the integrand of poles costing zero excitation energy, whose precise contribution depends 
on complicated geometrical factors which we have simplified. However the oscillatory 
r-dependence, the main focus of this work, remains unchanged from the free field result 


(12). This form is dictated by the final momentum integral always having a sharp lower 


limit p = p, arising from the physical situation of gapless excitations at a Fermi surface. 


Note also that the first correction (20) can potentially exceed the free-held result (12) for 


large px o, highlighting for the need for caution when interpreting perturbative estimates. 

In the next section the dependence of the wavefunction on a diquark source strength 
j will be investigated numerically using free lattice fermions. The term proportional 
to j in the action explicitly breaks baryon-number symmetry, and encourages diquark 
pairing which has the effect of opening up a gap at the Fermi surface; we will see that 
the impact on the oscillations is profound. 


3 Lattice Formulation 

The lattice results in this study employ two flavors of Wilson fermion, with action for 
the quarks given by pf] 

s = ££ 'lPixM X y(l-l)'IPiy KjS X y\ (C ^ T ^2x 75 ) ) (^l) 

x,y i=l,2 

with (in units where lattice spacing a = 1) 

M xy (p) =8 xv -kY j [(1 - 7 u)e^°U,(x)5y iX+ o + (1 + lv )e-^Ul{y)8 yjX .,\ . (22) 

V 

The subscripts on the quark fields denote flavor, and the Pauli matrix r 2 acts on color 
indices. The usual action is supplemented by the addition of gauge invariant scalar 
isoscalar diquark source terms with strengths j , here chosen real (gauge invariance fol¬ 
lows from the SU(2) identity T 2 UT 2 = U*). These terms serve a dual role; they enable the 
extraction of diquark condensates and anomalous quark propagators of the generic form 
('i/j(x)'iJj T (y)) , ('0 T (x)'0(i/)), and their inclusion also mitigates the long-wavelength fluc¬ 
tuations due to Goldstone modes in the superfluid phase which forms for p > p a — 
in which baryon number symmetry is spontaneously broken. As a direct consequence 
the number of iterations required to invert the matrix M is reduced by the opening of 
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Name 

Operator 

r 

I 

J p 

D°(x, r ) 


1 

0 

0+ 

D° 6 (x,r) 

ipi (x + r)KF ip 2 (x) 
—ip2(x + r)KT ipi(x) 

75 

0 

0- 

D^(x,r) 

li 

0 

1- 

7r(x, r) 


75 

1 

0- 

p{x,r) 


7* 

1 

1- 

Mp(x, r) 

ipi(x + r ) T ip 2 (x) 

1 

1 

0+ 



7s7 i 

1 

1+ 


Table 1: A summary of the diquark and meson operators used in the analysis, with K = C 75 T 2 . 


a gap in the Dirac spectrum of O(j ). In order to make contact with a theory where 
baryon number is a conserved charge we need ultimately to examine the limit j —> 0. 
The action (21) may be rewritten 


S = (ip,<f>) 


f M(n) 

\-Ktf7. 5 


*j 75 \ AA 

m(-h)J \<i>) 


VMifjiJ)*, 


(23) 


using the relabellings ip = ip 1 , ip = ipi, (p = —Cr 2 ipJ, (p = —ipJCT 2 - The enhanced matrix 
M obeys an important identity 

= (24) 


where T 5 = diag( 75 , 7 s). We can now outline hadron interpolation operators, following 
the treatment in [8j. In this study we will consider isovector meson excitations of the 
form 

Ml = ip\Yip2 = iPT(Ct 2 )P t (25) 

and baryonic diquarks in both isoscalar 

D° r = ipl (C75T 2 )T^2 = ip T ^V T 4> T = -0T 75 iP (26) 


and isovector 

D l v = ip T (C lb T2)ViP (27) 

channels. The Dirac matrices T and corresponding spacetime quantum numbers J p are 
listed in Tabled! 

The general quark propagator is 

<**) = M-'inJ) = S ^) . (28) 
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The relations (23p4) imply two distinct components: 


Normal: 

Anomalous: 


j) = M 

j ) ^(fnp 


\n,j) = Sj 

-/b -j)- 


5 {-V>-3)\ 


(29) 

(30) 


In the limit j —> 0 and on a finite lattice the anomalous components vanish. We 
can now write expressions for the hadron correlators required for the wavefunctions 
corresponding to spatial separation r between the (anti)quarks at the sink. First define 
eg. Mr(x,r) = ij)\{x)Yij} 2 {x + f), and = —^{x + r)T r ijji( y x) (with f = yoT^o). Then 

c mr( x ,y,r) = {Mr(y,$) M r(x,r)) = - y)T^S^(x + f, y^T } 

+ ^{S^X x + r,y)T^ 5 S^(y,x)^ 5 t}. (31) 

There are two terms, one involving a connected contraction of normal propagators, and 
one a connected contraction of anomalous propagators. 

Similarly in the isoscalar diquark channel (with D = (C757 2 )r , 0 2 and D 1 = 

'0 2 f(C75T 2 )^f): 

C&r ( x ,y,r) = { D r(y,®)) D r (x,r)) = + tv{S^(y,y)T^ 5 }tv{S^(x, x + r)T^ 5 } 

- tT{S^{x,y)TS^{x + f,y)T}. (32) 

Here the anomalous propagators only contribute to disconnected pieces, which we hence¬ 


forth ignore. A judicious use of (|24|) means the required normal contributions can be 

It is 


evaluated with a single inversion of Ai using a local source at y, just as in 
shown in [H] that the connected contribution to C° r is non-vanishing for T G { 1 , 75 , 7 ;}. 


Comparison of (31) and (32) also reveals the origin, in the limit y = j = 0, of the 
degeneracy between mesons with spin-parity J p and diquarks with noted in [ 8 ]. 
Once the system enters a superfluid phase for y > fi 0 baryon number is no longer a good 
quantum number, and we therefore expect degeneracy between meson and diquark chan¬ 
nels sharing the same J p . This prediction was approximately observed for the lightest 
states in [ 8 ], and must arise via non-trivial contributions from anomalous propagators. 
Finally, we note that a similar analysis for the isovector diquark correlator C \ r yields 


the same connected contribution as in (32), but this time with V restricted to ^ 757 ;. 


Next we examine wavefunctions evaluated for non-interacting quarks. Fig. [2] plots 
the wavefunction T}(F, r) defined by the ratio 


T(r,r) = 


^-C(T,r; 0 ;r) 

^-C(f,r;0;d)’ 


(33) 


where as usual the sum over the sink location x projects onto the zero momentum 
sector. The calculation is performed by inverting Ai on a 24 4 lattice, with n — | chosen 
to yield massless quarks. For consistency’s sake all curves shown correspond to temporal 
separation r = 8 between source and sink, and only points evaluated along a lattice axis 
are plotted to minimise the impact of discretisation artifacts. 
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Figure 2: Meson and diquark wavefunctions ’F(r, r = 8) for free massless quarks and various /r and j. 


Two trends are discernible. First, as chemical potential increases from 0 to na — 0.75 
the width of the main central feature of the free wavefunction decreases, and it develops 
an oscillatory structure whose spatial frequency is in qualitative agreement (but slightly 


smaller) than the prediction of (12). The amplitude of the oscillations dies away more 
quickly than for those observed in the 2 + Id Gross-Neveu model [15] , since in 3 + Id the 
profile is given by Ji(r/r) rather than J 0 (rfi). Second, as diquark source j is increased 
at fixed /i, both the width of the central feature (especially for /i = 0) and, for /j, > 0 the 
amplitude of the oscillations are diminished. This can be understood as a consequence 
of the opening of a gap at the Fermi surface over a momentum range A O(j), as 
shown in Fig. [lj with the consequence that there is mixing between particle and hole 
states so that both higher momentum states contribute to the correlators, and the sharp 


momentum cutoff in the integral leading to (12) is smeared out 


4 Numerical Results 


Hadron wavefunctions formed from interacting quarks were calculated using QC 2 D en¬ 
sembles generated using the quark action ( 2lf22 ) together with an unimproved Wilson 
gauge action for the gluons. The simulation parameters were /3 = 1.9, k = 0.168, 
corresponding to lattice spacing a = 0.178(5)fm with scale set by assuming the string 
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tension is (440MeV) 2 , and mass ratio m n /m p = 0.807(5) [5|. Most results are obtained 
on a 12 3 x 24 lattice, corresponding to a physical temperature T = 44(2)MeV, although 
for pa = 0.25 we also have results from 16 3 x 24 for comparison. This temperature is 
sufficiently low to support the existence of an extended range of p in which the the¬ 
ory is simultaneously confining (as indicated by a near-vanishing Polyakov loop) and 
superfluid (as indicated by a non-vanishing condensate (ipj (£ 75 ) 72 ^ 1 ) 7 ^ 0 as j —» 0 ). 
With these parameters values of pa in the range [0.0,1.1] were explored; the onset value 
p D = = 0.323(3)a _1 . The so-called “quarkyonic” regime where baryon density, 

pressure and superfluid condensate all scale with p according to the expectations of a 
system of degenerate quarks, lies approximately in the range pa G (0.4, 0.8) [5', [ 6 ]. 

For |r| > 0 the point-split correlators ( 3lp2 ) are not gauge invariant without an 
insertion of path-ordered link variables along some selection of paths joining the two 
halves of the sink. To mitigate the effects of the signal fluctuations introduced by this 
non-unique procedure, we instead choose to gauge-fix the configuration and use unit 
links to complete the loop. We fix a discretised Coulomb gauge defined by 


A °w = Xi bfw - A ?( x -*)] = 0, 


(34) 


i =1 


where the gauge transformation G(x) extremises the functional 


F[U a ] = -ReTY^£(7fl 


X 


(35) 


i=l 


with U^(x) = G(x)U p (x)G~ 1 (x + p). To achieve this, the simplest algorithm [TF] one 
can adopt is a local procedure which visits one lattice site at a time and attempts to 
minimize its contribution to the functional ( |35j ), which can be written as: 

F\ oc (x) oc -Re Tr ^ [U p (x) + U p (x - p)] . (36) 


Two observables are usually monitored during this procedure. One is the functional (35) 
itself, which decreases monotonically and eventually reaches a plateau. The other one 
is a measure of the first derivative of F[U] during the gauge-fixing process defined as 


0 a =4 JjTr(A G a)(A G )'a)]. 

X 


(37) 


where V is the lattice volume. This quantity eventually approaches zero when F[U] 
reaches its minimum and can be used as a stopping parameter for the procedure. Here 
we chose 6 < 10 _3 °. 

As shown in ra. the gauge-fixing procedure is affected by long-range correlations 
when the size of the lattice increases, which eventually leads to a critical slowing down. 
There are a number of strategies available to improve this and here we adopt the so-called 
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Figure 3: Meson (left) and diquark (right) Euclidean correlators evaluated for various /u at ja = 0.02. 


overrelaxation procedure [TS]. This method replaces the gauge fixing transformation 
G(x) with its power G u (x ) at each iteration, where u G (1,2) is a parameter which is 
tuned empirically at an optimal value, depending on the volume of the lattice and the 
coupling /3. The exponentiation is obtained by a truncated binomial expansion 


N 


G“ = £ 


7n(Wj 


n =0 


n\ 


(G — I) n , 7 nH 


T(u + 1) 

r (u + 1 — n) ’ 


(38) 


with 2 < N < 4. ffere we chose u ; = 1.75 and N — 4. Before G“{x) is applied to the 
link, it has to be reunitarised in order to belong to the gauge group. 

Fig. [ 3 ] shows typical meson and diquark timeslice correlators ^C mi i(0,a:;r = 0) 
evaluated for /ia G [0,1.0). By construction mesons yield a signal which is symmetric 
about the temporal midpoint of the lattice, whereas diquarks carrying baryon charge 
are time-asymmetric once /1 > 0. We interpret the forwards-moving lighter states as 
diquarks and the backwards-moving heavier states as anti-diquarks (in contradistinction 
to the nomenclature adopted in i)- Inspection of the plots reveals that while the 
1 + meson is more massive than the 1“ diquark at // = 0 (recall they should only be 
degenerate in the limit j —>• 0), both become lighter as /j increases beyond onset. 

This is in accordance with the spectrum results found in [8] using all-to-all propaga¬ 
tors evaluated using source dilution. In brief, in mesonic channels, the 0“ pion and 1“ 
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rho are the lightest states at fj, = 0, but beyond onset there is a level crossing, the rho 
becoming lighter while the pion mass increases almost linearly with . The two lightest 
states for fi > fi Q are, however, the isoscalar 0 + and the slightly heavier isovector 1 + . 
In the diquark channel, the isoscalar 0 + and isovector 1 + are approximately degenerate 
with the pion and rho states at /i = 0 but, again, beyond onset become much lighter and 
approximately degenerate with the mesonic states sharing the same quantum numbers, 
as expected in a superfluid ground state. In the same post-onset regime isoscalar 0“ and 
1“ states are observed at roughly the same mass scales as their tt/ p counterparts. 

The main focus of this paper is the post-onset regime /i > /i G , where we expect 
degeneracy between meson and baryon states due to the superfluid nature of the ground 
state, so that physical states are labeled solely by their J p quantum numbers. However, 
meson and baryon operators need not have the same overlap with the actual bound 
states; Table 2 of Ref. [8] shows the quality of the signal for each operator. We will 
mainly concentrate on the channels yielding a “good” signal for p, > fi 0 : the 0 + , 0“ and 
1“ diquarks and the 1 + meson. 

By r = 8 the temporal decay in both cases is well-approximated over a wide range 


of r by a single exponential, indicative of a bound state. In this regime the ratio (33) 
can be written eg. 


T(r, r) = 


E» I]x(°l^(°) r ^(o)N)(^( 3: ) r ^( :r + me EiX0 _ +rQjo) 

E 1 Ex(°IV i ( 0 ) r V'(0)|i)(^(a:)r^(a:)|0)e 


-Eix o 


(39) 

where | H) is understood as the lightest hadron state in the particular channel, and can 
be seen to be proportional to the wavefunction defined in eqn. (|2]) but normalised to 
unity at the origin r — 0. Unless otherwise stated, we will show results for wavefunctions 
T evaluated at this temporal separation. 

Fig.g shows wavefunctions T (r, r) evaluated on timeslice 8 for the same four chan¬ 
nels of Fig [ 2 } The plots compare data from 12 3 x 24 and 16 3 x 24 lattices; data from 
all available separations r up to half the lattice extent is shown. There is some evi¬ 
dence for discretisation artifacts, eg. around r ~ 2a, but overall the trend with r is 
smooth within errors; moreover finite volume effects are absent for spin-0 states, and 
larger (but still consistent with zero) for spin-1, where data from the smaller volume are 
considerably noisier. The absence of spatial volume effects is consistent with excitations 
being confined states, with little or no contamination from image charges which might 


be anticipated from the algebraic fall off (12) of the free-held wavefunctions. 


Fig- 1] shows wavefunctions calculated with varying r values. For clarity’s sake only 
data from separations r on-axis are plotted. The 1“ diquark and the 1 + meson show 
similar behaviour; the approximate consistency of data from different r is good evidence 
for the existence of a discrete bound state and the applicability of the approximation 


(39), although for /ra < 0.5 the data is much nosier for the meson than the diquark. 


For n < Ho the wavefunction has a large width and has not yet vanished by the lattice 
midpoint, which is consistent with the volume dependence of the spin-1 data seen in 
Fig. |U 
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Figure 4: Hadron timeslice correlators evaluated for ja 


0.02 and fia = 0.25 for various channels. 


In Fig. [6] we show the consequences for the scalar diquark wavefunction of varying 
the diquark source j\ the data have been generated in a partially-quenched approach 
using an ensemble generated with ja = 0.02. It can be seen that for fi > 0 increasing 
j has the effect of very slightly shrinking the wavefunction. In contrast to the data for 
free fermions in Fig[2j there is no sign of oscillatory behaviour developing in the j —> 0 
limit. 

Figs. [7] show r = 8 wavefunctions for the four channels of interest as // is varied. 
Comparison of the plots enables us to order the states by spatial size at /i = 0: 0 + < 1“ < 
1 + . The 0“ data are considerably noisier. The main common trend is the systematic 
decrease in the size of all states as /j increases, so that by /ia ~ 0(1) the maximum extent 
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Figure 5: Comparison between diquark (top) and meson (bottom) wavefunctions evaluated at different 
timeslices values r = 7,8,9,10. We show the results for 4 different values for the chemical potential 

At = 0.0,0.25,0.5,0.75. 
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Figure 6: 


Diquark wavefunctions relative to Dj showing the j-dependence for various values of the 

chemical potential jj,. 


ra ~ 3. With the exception of the noisy 0 _ , there is no sign of any of the wavefunctions 
changing sign or developing oscillatory behaviour as /i increases. 

To crudely quantify the evolving spatial size, we extracted from the data the full 
width at half the maximum of T(r, r), which was obtained by fitting the wavefunction 
to a spline. In Fig. [8] we then plot the resulting width c(/u) in each channel. The 
results are plotted both as a function of fi and of l//q and confirm the trends reported 
above, and also highlights that post-onset cr(0”) is larger than both 0 + and 1 + states. 
We note that cr(0 + ) > <r(l _ ) for fia > 0.6, which is confirmed by close inspection 
of Fig. [7J For fia > 0.6, or 1/fia < 1.5, the widths approach one another making the 
various channels difficult to distinguish, but the right panel of Fig. [8]suggests that bound 
states are increasingly dominated by a single length scale a oc // -1 , with the hierarchy 
cr(l“) < O r (0 + ) < cr(l + ) < cr(0 - ). 

5 Discussion 

We have presented results from the first attempt, using orthodox lattice techniques, to 
examine the spatial structure of gauge-invariant excitations in a baryonic medium. The 
results complement a previous study [8] of the excitation spectrum. We have focused 


18 


















































































^(r, r) 



r 


Figure 7: Diquark and meson wavefunctions for various values of the chemical potential /i at ja = 0.02. 
Below onset, some operators are very noisy and are excluded from the plots for clarity. 
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Figure 8: Estimated width a of the wavefunction for the different operators, obtained by fitting to a 
spline and then extracting the width at half maximum, as a function of y (left) and yT 1 (right). 


on the post-onset regime /i > /x 0 in which baryon charge density is non-zero in the 
T — y 0 limit, and the ground state is a superflnid. Onr results are consistent with 
the indistinguishability of mesons and diquarks in a superfluid, and suggest a scale 
hierarchy cr(0 + ) ~ cr(l _ ) < cr(0~) < a(l + ), to be compared with the mass hierarchy 
m(0 + ) < m( 1 + ) <C m(l _ ) < m(0 _ ) found in [8j. As a general rule, signals obtained in 
diquark channels were less noisy than those from mesons. 

The scale hierarchy becomes less well-defined as /x increases and the wavefunctions 
shrink; from /xa ~ 0.6 onwards all the channels yield wavefunctions of approximately 
equal extent. For orientation, if the string tension is used to set the scale this corresponds 
to fj, ~ 670MeV, at which point the quark density n q ~ usb ~ 5fm -3 (using Q), or 
roughly 10 x nuclear density [5]. Spline hts to the profiles then yield the approximate 
behaviour a oc /x -1 with a different hierarchy a(l _ ) < cr(0 + ) < cr(l + ) < cr(0 _ ). This is 
consistent with the expectation cr oc /x” 1 which assumes that /x ~ kp is the only relevant 
scale at high density. Indeed, this is precisely the content of the free-held prediction (12). 
The physical picture is that bound-state excitations are formed from quarks close to the 
Fermi surface with a characteristic de Broglie wavelength Xp ~ jtx -1 . The absence of 
appreciable finite volume effects suggests, however, that the influence of image charges is 
negligible, and that confinement continues to hold. The conjunction of both properties 
characterises the so-called quarkyonic regime. 

Finally, the absence of oscillatory behaviour in T(r, r) at low temperature, in contrast 
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with the weak-coupling prediction ( 12 ) needs some consideration. One obvious departure 
from weak-coupling behaviour is the formation of a superfluid condensate (qq) 7 ^ 0 , 
which for a degenerate system should, via the BCS mechanism, induce an energy gap 
A ~ Aqc 2 d ~ {qq) /7-t 2 at the Fermi surface. The presence of a gap removes the sharp 
momentum cutoff in the integrals leading to the expressions ( 12 , 20 ), and therefore also 
the oscillations. A gap can also be modelled in free-held theory by the introduction of a 
diquark source j 7 ^ 0, and the curves shown in Fig. [2] strongly suggest this explicit gap 
does indeed dampen the oscillations. The fact that oscillations remain absent from the 
wavefunctions of interacting quarks as j —y 0, demonstrated in Fig. [ 6 j is consistent with 
the post-onset QC 2 D gap being generated dynamically. 
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